#include <iostream>
#include "includes/slope.hh"
double Slope::minmod(double a, double b){
  double tmp=0.0;
  if (a>=0 && b>=0)
    tmp=std::min(a,b);
  if (a<=0 && b<=0)
    tmp=std::max(a,b);
  return tmp;
}

double Slope::minmod3(double a, double b, double c){
  double tmp=0.0;
  if (a>=0 && b>=0 && c>=0)
    tmp=std::min(std::min(a,b),c);
  if (a<=0 && b<=0 && c<=0)
    tmp=std::max(std::max(a,b),c);
  return tmp;
}

double Slope::maxmod(double a, double b){
  double tmp=0.0;
  if (a>=0 && b>=0)
    tmp=std::max(a,b);
  if (a<=0 && b<=0)
    tmp=std::min(a,b);
  return tmp;
}

Slope::Slope(){
  slopeleft.erase(slopeleft.begin(),slopeleft.end());
  sloperight.erase(sloperight.begin(),sloperight.end());
  slopecenter.erase(slopecenter.begin(),slopecenter.end());
  limitslope.erase(limitslope.begin(),limitslope.end());
  varleft.erase(varleft.begin(),varleft.end());
  varright.erase(varright.begin(),varright.end());    
  std::vector<std::vector<double> >(slopeleft).swap(slopeleft);
  std::vector<std::vector<double> >(sloperight).swap(sloperight);
  std::vector<std::vector<double> >(slopecenter).swap(slopecenter);
  std::vector<std::vector<double> >(limitslope).swap(limitslope);
  std::vector<std::vector<double> >(varleft).swap(varleft);
  std::vector<std::vector<double> >(varright).swap(varright);  
  slopeleft=std::vector<std::vector<double> >(memory::phynumber, std::vector<double>(gridsetting::dimension));
  sloperight=std::vector<std::vector<double> >(memory::phynumber, std::vector<double>(gridsetting::dimension));
  slopecenter=std::vector<std::vector<double> >(memory::phynumber, std::vector<double>(gridsetting::dimension));
  limitslope=std::vector<std::vector<double> >(memory::phynumber, std::vector<double>(gridsetting::dimension));
  varleft=std::vector<std::vector<double> >(memory::phynumber, std::vector<double>(gridsetting::dimension));
  varright=std::vector<std::vector<double> >(memory::phynumber, std::vector<double>(gridsetting::dimension));
  for(int tmpp=0; tmpp<memory::phynumber; tmpp++)
    for(int tmpd=0; tmpd<gridsetting::dimension; tmpd++){
      slopeleft[tmpp][tmpd]=0.0;
      sloperight[tmpp][tmpd]=0.0;
      slopecenter[tmpp][tmpd]=0.0;
      limitslope[tmpp][tmpd]=0.0;
      varleft[tmpp][tmpd]=0.0;
      varright[tmpp][tmpd]=0.0;                  
    }
};

Slope::~Slope(){
};

void Slope::get_slope(SubprojectBase* subproject,double* globalcoord, double* unitouternormal,double pillarlen,int i,int j, int currentdim,int* interscounter){
  int interscountersum=0;
  double slopetmp1,slopetmp2;
  double vartmp[memory::phynumber];
  double dummyflag=0.0;
  double var_i_middle[memory::phynumber];
  double var_j_middle[memory::phynumber];
  basis.getVarMiddle(i,var_i_middle);
  
  
  for(int tmpi=0;tmpi<gridsetting::dimension;tmpi++)
    interscountersum+=interscounter[tmpi];  

  // internal left/right_variable
  if(j>=0){
    basis.getVarMiddle(j,var_j_middle);    
    for(int tmpi=0; tmpi<memory::phynumber; tmpi++){
      if(unitouternormal[currentdim]<0){
	varleft[tmpi][currentdim]=var_j_middle[tmpi];
      } else {
	varright[tmpi][currentdim]=var_j_middle[tmpi];
      }
    }
  }
  // bondary left/right_varible
  if(j<0){
    subproject->boundary->boundary_function(globalcoord,i,unitouternormal[currentdim],vartmp,dummyflag);
    for(int tmpi=0; tmpi<memory::phynumber; tmpi++){
      if(unitouternormal[currentdim]<0){
  	varleft[tmpi][currentdim]=vartmp[tmpi];
      } else {
  	varright[tmpi][currentdim]=vartmp[tmpi];	
      }
    }
  }
  // left/right_slope  
  for(int tmpi=0; tmpi<memory::phynumber; tmpi++){
    if(unitouternormal[currentdim]<0){
      slopeleft[tmpi][currentdim]=(var_i_middle[tmpi]-varleft[tmpi][currentdim])/pillarlen;
    } else {
      sloperight[tmpi][currentdim]=(varright[tmpi][currentdim]-var_i_middle[tmpi])/pillarlen;
    }
  }
  // center_slope
  if(interscounter[currentdim]==2)
    for(int tmpi=0; tmpi<memory::phynumber; tmpi++){
      slopecenter[tmpi][currentdim]=(varright[tmpi][currentdim]-varleft[tmpi][currentdim])/(2.0*pillarlen);
    }
  // limited_slope
  if(interscountersum==2*gridsetting::dimension){
    for(int tmpi=0; tmpi<memory::phynumber; tmpi++){
      for(int dimi=0; dimi<gridsetting::dimension; dimi++){
        switch(setting::limiter) {
    	case setting::minmod:
     	  limitslope[tmpi][dimi]=minmod(slopeleft[tmpi][dimi],sloperight[tmpi][dimi]);
        break;
        case setting::superbee:
    	  slopetmp1=minmod(sloperight[tmpi][dimi],2.0*slopeleft[tmpi][dimi]);
    	  slopetmp2=minmod(2.0*sloperight[tmpi][dimi],slopeleft[tmpi][dimi]);
          limitslope[tmpi][dimi]=maxmod(slopetmp1,slopetmp2);	  
        break;
        case setting::mc:
    	  limitslope[tmpi][dimi]=minmod3(slopecenter[tmpi][dimi],2.0*sloperight[tmpi][dimi],2.0*slopeleft[tmpi][dimi]);
        break;
	}
      }
    }
    // update slope
    //    memory::variable[memory::phynumber][i]=slopeleft[0][0]*pillarlen;
    //memory::variable[memory::phynumber+1][i]=slopeleft[1][0]*pillarlen;
    // memory::variable[memory::phynumber][i]=slopeleft[0][0]*pillarlen;
    // memory::variable[memory::phynumber+1][i]=sloperight[0][0]*pillarlen;
    // memory::variable[memory::phynumber+2][i]=slopecenter[0][0]*pillarlen;
    // memory::variable[memory::phynumber][i]=limitslope[0][0]*pillarlen;
    // memory::variable[memory::phynumber+1][i]=limitslope[1][0]*pillarlen;        
    
    for(int tmpi=0; tmpi<memory::phynumber; tmpi++){
      for(int dimi=0; dimi<gridsetting::dimension; dimi++){
    	memory::limit_slope[tmpi][dimi][i]=limitslope[tmpi][dimi];
	// memory::variable[memory::phynumber+tmpi*dimi][i]=memory::limit_slope[tmpi][dimi][i];
	// memory::variable[memory::phynumber+tmpi][i]=sloperight[tmpi][dimi]*pillarlen;
      }
    }
  }
};
